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Abstract. We consider the quantum simulation of relativistic quantum 
mechanics, as described by the Dirac equation and classical potentials, in trapped- 
ion systems. We concentrate on three problems of growing complexity. First, we 
study the bidimcnsional relativistic scattering of single Dirac particles by a linear 
potential. Furthermore, we explore the case of a Dirac particle in a magnetic 
field and its topological properties. Finally, we analyze the problem of two Dirac 
particles that are coupled by a controllable and confining potential. The latter 
interaction may be useful to study important phenomena as the confinement and 
asymptotic freedom of quarks. 
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1. Introduction 

The field of quantum simulators is one of the most rapidly growing in quantum 
information science. It was Feynman [T] who stated thirty years ago that a controllable 
quantum device could emulate the dynamics of another quantum system exponentially 
faster than a classical computer. Since then, this hypothesis has been confirmed [2] 
and, subsequently, intensive theoretical and experimental research has followed [2l|4]. 

Quantum simulators of quantum relativistic systems are among the most elegant 
and, among a few others, the most promising candidates for going beyond classical 
computational capabilities. In the last years, the simulation of black holes in Bose- 
Einstein condensates (BEC) [5] and properties of the expanding universe [OJ [7J have 
been proposed, and an experiment on BEC sonic black hole has been realized [Hj. A 
relevant milestone was given by the proposal for the simulation of Dirac equation and 
associated quantum relativistic effects in a single trapped ion [9] , that was subsequently 
experimentally realized |10j . This was the first experimental observation of the physics 
associated with the Zitterbewegung phenomenon, the fast quivering motion stemming 
from the Dirac equation and predicted by Schrodinger in the early days of quantum 
mechanics [11]. Further developments produced the theoretical proposal 12 and 
experimental realization [13] of the simulation of Klein paradox 14, and a proposal for 
simulating the Majorana equation and unphysical operations like charge conjugation 
or time reversal with trapped ions [15] . 

Many other proposals and some experiments, in a wide variety of systems, 
have recently appeared, as the simulation of Zitterbewegung in semiconductor 
quantum wells [16] or in graphene [17l [18], Klein paradox in graphene [19], Dirac 
oscillator in a trapped ion [20], Zitterbewegung and Dirac physics with ultracold 
atoms [2T1 (22J [231 1211 (2SJ [20J, Klein paradox with atomic ensembles [27], optical 
Zitterbewegung in metamaterials [28U29) . derealization of relativistic Dirac particles in 
cold atoms |30) . photon wave function and Zitterbewegung [31 , similarity of electron's 
Zitterbewegung to the Adler-Bell-Jackiw anomaly in QED and its manifestation in 
graphene [32], photonic analog of Zitterbewegung in binary waveguide arrays [33j . 
Zitterbewegung theory in multiband Hamiltonians 34J, classical Zitterbewegung 
in reduced plasma dynamics [35) . Zitterbewegung analogs in nonlinear frequency 
conversion [36j , experimental realization of an optical analog for relativistic quantum 
mechanics in an optical superlattice (37) . relation between parity and Zitterbewegung 
and proposed simulation in trapped ions |38j , Zitterbewegung in a magnetic field and 
proposal for trapped ion simulation [39) , Wilson fermions and axion electrodynamics in 
optical lattices [JO] , the Schwinger effect for a possible implementation with atoms in 
optical lattices [JJJ , Dirac equation for cold atoms in artificial curved spacetimes [32] , 
a theoretical analysis of cold atom simulation of interacting relativistic quantum field 
theories [43] , or an analysis of the photonic simulation of the quark model [44] . For a 
review of Zitterbewegung of electrons in semiconductors, see Ref. |45) . 

In this article, we make a step forward and analyze novel features of relativistic 
quantum mechanics simulations with trapped ions. We first analyze Klein paradox in 
2+1 dimensions for different kinds of potentials. We address the problem of 2+1 Klein 
tunneling which contains novel and interesting features, like entanglement between 
the transmitted/reflected wave packets in the Klein dynamics and the transverse 
momentum wave function. We also point to the fact that more general nonlinear 
V(x, y) potentials will already demand computational classical resources that are 
approaching what is currently feasible. Thus, a general 2+1 Klein dynamics will 
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already be an interesting problem to be addressed by a quantum simulator. Then, 
we analyze a Dirac particle in a magnetic field, making here an emphasis on its 
topological properties, and show how it can be simulated with trapped ions. Finally, we 
study the quantum simulation of two Dirac particles coupled by a confining potential 
implemented in a system of three trapped ions. This model represents a simplified 
version of the MIT bag model of nuclear physics [46l [47] . This semianalytic approach 
is still frequently used for analyzing quantum chromodynamics (QCD) in the non- 
perturbative regime. In this manner, we show that one could implement interesting 
QCD features like asymptotic freedom and confinement in a table-top experiment, 
controlled just by turning on and off a laser. 

The paper is organized as follows. In Section [2j we briefly review the quantum 
simulation of the Dirac equation in trapped ions. In Section [3j we analyze the Klein 
paradox in 2+1 dimensions. In Section [4] we study the simulation of a Dirac particle 
in a magnetic potential and its topological properties. In Section [5] we address the 
problem of two Dirac particles coupled via a confining potential simulated with three 
trapped ions. Finally, we present our concluding remarks in Section [6] 

2. Dirac equation simulation in trapped ions 

The Dirac equation is arguably the most fundamental wave equation [48] . It 
is historically considered as a significant step forward towards unifying quantum 
mechanics and special relativity, accurately describing the hydrogen atom spectrum, 
while predicting ab initio spin and antimatter. The Dirac equation acquires full 
significance in the context of quantum field theory and second quantization, where the 
number of quanta is not fixed. On the other hand, at the level of relativistic quantum 
mechanics, the Dirac single-particle solutions predict already intriguing phenomena. 
Among them, the best known are the Zitterbewegung |11| and Klein paradox (14) . 
The 3+1 Dirac equation reads 

ih^ =H^ rl ip = (ca-p + I3mc 2 )4>, (1) 

where a and /3 are the Dirac matrices that obey the Clifford algebra structure, 
{ai,ctj} = 2dij, {oti,P} = 0, /3 2 = I4, that appear when linearizing the expression 
for the relativistic energy, E — ^Jp 2 c 2 + m 2 c 4 . 

The Hamiltonian operator in Eq. ([I]) expressed in its "supersymmetric" 
representation [48], takes the form 

1/3+1 _ ( c(ff ■ p) - imc 2 \ , 

rt ° ~ ^ c(a ■ p) + imc 2 J ' [ > 

where a := (a ) = off-diag((7, (?) is the velocity operator, a := ((T x ,a y ,(7 z ) are 

the Pauli matrices, f3 := off-diag(— H2, H2), c is the speed of light and mc 2 the electron 
rest energy. 

A quantum simulation of the Dirac equation is a useful playground to analyze 
different regimes and to verify in a table-top experiment predicted phenomena like 
Zitterbewegung [9] HQ] or Klein paradox [12l [13] . To perform just the Dirac equation 
and Zitterbewegung, one considers a single trapped ion of mass M inside a Paul 
trap, with motional-mode frequencies v x , v y , and v z [S]. In the 3+1 case, the four- 
component Dirac bispinor will be codified in four internal levels of the trapped ion, 
\a), \b), |c), and \d). 
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For any spacetime dimensions, one should realize that the Dirac equation contains 
basically couplings cr^ (spin-orbit), and mc 2 aj (mass term). In general, one may 
implement the first kind of couplings, o~iPi, by a combination of Jaynes-Cummings 
(JC, also named red-sideband) and anti-Jaynes-Cummings (AJC, also named blue- 
sideband) interactions. The JC interaction excites one quantum of vibration while 
deexciting the internal state of the ion. The AJC interaction, in turn, excites one 
quantum of vibration while exciting the internal state of the ion. These interactions 
may be obtained by suitably chosen lasers, either in Raman or quadrupole-transition 
configurations, depending on the chosen ion. 

The resonant JC Hamiltonian can be written as 

H T = HrjCHa+ae*** + ^aV^), (3) 

where Q is the coupling strength, a + and a~ are the raising and lowering spin-1/2 
operators, and a and a) are the annihilation and creation operators associated with 
the corresponding motional degree of freedom, either x, y, or z. r/ — ky/h/2Mv is 
the Lamb-Dicke parameter [49] , where k is the wave number of the driving field. The 
AJC Hamiltonian reads 

H h = hi^a+a^e 1 ^ + a~ ae'^). (4) 

By properly adjusting laser phases and frequencies, one may combine a JC and 
an AJC interactions to obtain the kind of couplings aimed for, namely 

H%* = ihrj x Q x a x (al - a x ) = 2i] x A x n x a x p x , (5) 



with i(a x — a x )/2 = A x p x /h, where A x := ^Jhj2Mv x is the spread in position 
along the x-axis of the ground state harmonic-oscillator wavefunction and p x the 
corresponding dimensioned momentum operator. These interactions may be addressed 
to different internal levels and motional modes, in order to get the terms in Eq. ([2]) 
that are linear in the momenta. 

For implementing the second kind of couplings, mc 2 Oj, one may consider the 
carrier interaction between two internal levels of the ion. This consists of a coherent 
excitation of the internal level, while leaving the motion unaffected. The carrier 
Hamiltonian is 

H a = hflia+e^ + (J-e"^), (6) 

that, for appropriately chosen phases, gives rise to the term Mla y: needed for the mass 
terms for 3+1 dimensions. 

For the 3+1 dimension simulation [5], notice that the Hamiltonian 

Hl +1 = 2 V An(a a x d + a b x c )p x + 2rjAft« - a b y c )p y 

+ 2 V An(a a x c - a b x d ) Pz + Ml« + of). (7) 

is composed of basic JC, AJC and carrier units as exposed above. When this is 
expressed in matrix form, in the basis \a), \ b), |c), and \d), 



R 3+i_f 2riA£l(a-p)-ihQ, 



(8) 



it coincides with the Dirac equation Hamiltonian in 3+1 dimensions, Eq. with 
the equivalences for the speed of light and rest energy, 

c := , mc 2 := Ml (9) 
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An important property of quantum simulations of relativistic quantum mechanics 
with trapped ions is that the values of the speed of light and rest energy in Eqs. (9) 
may be controlled at will just by changing the laser intensities, f2 and fl. Thus, one 
may explore the transition from massless to massive Dirac fermions, or, equivalently, 
the transition from ultrarelativistic to nonrelativistic physics. 

One of the most astonishing predictions of the single free-particle solutions of the 
Dirac equation is the fast quivering motion called Zitterbewegung. It is unexpected 
because it predicts an oscillatory motion of a freely propagating electron. Thus, 
Galileo's inertia law is not fully verified for a free relativistic electron, contrary to 
the Schrodinger equation case, and the free Dirac particle is expected to quiver 
around in the absence of potentials. The reason for this is the non-commutativity 
of the components of the velocity operator, which is given by ca. Thus, the marriage 
between quantum mechanics and special relativity seems to contradict the inertia law, 
at least to some extent and within the subtle realm of relativistic quantum physics. 
The Zitterbewegung phenomenon has not been observed so far for a real relativistic 
electron, given that the predicted frequency, ~ 10 21 s _1 , and amplitude, ~ 10~ n cm, 
are difficult to access experimentally. Moreover, its correct physical prediction and its 
validity has been constantly questioned in the last decades, whether we remain in the 
domain of relativistic quantum mechanics or quantum field theory. 

The expression for the Dirac electron's position operator r — (x, y, z) in the 
Heisenberg picture, derived from df/dt = [r,Hn]/ifi, for 3+1 dimensions, reads [H] 

m = m + ^ , + f (a - j* ) (e*« + V» - l) , 

(10) 

where the first term on the r.h.s. is the initial position, the second is just the inertia 
law term, and the third one is the term associated to the Zitterbewegung. 

In the trapped ion simulation, the position of the ion mimics the position of 
the simulated Dirac particle. The ionic position operator evolution under the Dirac 
dynamics in Eq. (J8|, in 3+1 dimensions, is 

m = m + H*™£ t+ ( a - 2 0) - 1) . 

(ii) 

and the Zitterbewegung frequency can be estimated in the simulation as 



w ZB ps 2\E D \/H = 2^Ar 1 2 A 2 tt 2 p 2 /h 2 + il 2 « 2\J Nr] 2 n 2 + il 2 , (12) 

where E^> = (Hd) is the average energy, po is the average momentum for a peaked 
distribution, and ./V = (a^a) is the phonon number, respectively. Additionally, one 
can estimate the Zitterbewegung amplitude to be 



mc 



2 



rjh 2 nnA 



RzB 2mc V E ) A^VP^pl + h 2 ® 2 ' ^ 

and i?ze ~ A, if rjfl ~ f2, which can be measurable in an experiment, as was shown 
in the 1+1 case in Ref. [TO"] . 
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3. Bidimensional relativistic scattering simulated in trapped ions 

The Klein paradox is another counterintuitive prediction of the relativistic quantum 
mechanical solutions of Dirac equation. A Dirac particle may tunnel through a steep 
barrier and propagate indefinitely regardless of the fact that the potential may be 
extended until arbitrarily long distances. The reason for this is that the positive energy 
electron may become, upon collision with the barrier, a negative energy electron and 
as such propagate freely throughout the barrier, wherever the potential barrier is 
V > 2mc 2 . The standard explanation in quantum field theory is that there is enough 
energy in the system to produce a particle/antiparticle pair. 

The theoretical proposal [H] of the simulation of the Klein paradox in 1+1 
dimensions was based on the Hamiltonian 

rt D = ca x p x + mc a z + ax, (14) 

where a is the potential gradient constant. The experimental simulation of the Klein 
paradox relied on the manipulation of two trapped ions [T3] . The center of mass mode 
together with the internal degrees of freedom of one of them simulates the positive 
and negative energy spinors, while the second one implements the external potential 
that produces the Klein paradox [TH [T3] . In order to measure the motional state of 
the ions in the previous simulations, the standard approach is to couple the motion to 
the internal states of one ion and, then, apply fluorescence detection to measure the 
internal state with novel techniques (50) [51] [52j [Qj [10] . 

Here, we propose an extension of the Klein paradox simulation to 2+1 dimensions. 
We will show that, for some potentials, this problem may be addressed with the 
techniques already developed for 1+1 dimensions. On the other hand, for arbitrary 
potentials in 2+1 dimensions, one is already approaching the limit of classical 
computational power and, in this quantum simulator could already predict 

physics beyond current capabilities. At the end of this section, we give a heuristic 
analysis of this issue. 

We begin by analyzing the Dirac equation in 2+1 dimensions with a linear 
potential along the x coordinate and constant along the y coordinate. The expression 
reads 

H D — ca x p x + c<TyPy + mc a z + ax, (15) 

where a is the potential gradient constant. In order to analyze this dynamics, one can 
point out that p y is a constant of motion: [H^ l:2+1 ,p y ] = 0, such that one may solve 
the problem for each value of p y . One can realize that ca y p y + mc 2 a z , with p y fixed, 
can be expressed as a new Pauli matrix with an effective mass coefficient, 

cOyP y + mc 2 a z = rhc 2 o y , (16) 

with mc 2 = \jpyC 2 + m 2 c 4 . Here ay = n y a y + n z a z is an effective Pauli matrix in 

the y ~ z plane (with a basis change in that plane), with (n y ,n z ) — (p y /mc,m/m). 
Accordingly, the problem, for each p y , is reduced to a Klein Hamiltonian in 1+1 
dimensions, 

^D,cff +1 (Py) = C(J xPx + mc 2 a y + ax. (17) 

In Fig. [T] we illustrate the dependence of the Klein tunneling on the incoming 
momentum orthogonally to the potential gradient, p y . To do so, we simulate a 
wavepacket which initially has a Gaussian profile on the x direction and a well defined 
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momentum p y . We plot the resulting wavefunction (normalized just in x and constant 
in p y for the sake of comparison between different p y components), \ip(x,p y ,t)\ , for 
different time snapshots t = 0, 20, 40, 60, 80, 100, in units oi h = c — a = 1, m = 0.5. 
These plots evidence the rapid suppression of Klein tunneling for increasing p y . Given 
that the transmitted amplitude decreases exponentially with the effective mass m as 
exp(—nfh 2 c /hca) [12] . we find that already at p y = 1 almost all the wavefunction 
gets reflected. 

The plots in Fig. [T] admit another interpretation, based on the entanglement 
between the transmitted and reflected wave packets in the x coordinate and the 
transverse momentum along the y direction. In the language of quantum field theory, 
the electron-positron pair production along x direction will be conditional to the 
specific value of each momentum wave packet component along y direction. Thus, 
one could control pair production just by changing incidence angle, p y /p xi and even 
perform quantum logic on pair creation conditional to the value of the transverse 
momentum. 

Based on the behaviour at fixed values of p y we can reconstruct the dynamics of 
arbitrary wavepackets, fully in position space. More precisely, we can write 

m))=J2 I dp y dxe- im ^ +1 ( p yV h Mx,Py>s)\x)®\p y )®\s), (18) 

s J 

where tpa(x,p y , s) is the initial state of the wavepacket, expressed in position space 
for the X coordinates, \x), momentum for the Y coordinate, \p y ), and spinor state, 
\s = 0, 1). This reconstruction formula evidences the conditional evolution mentioned 
before, and in particular the entanglement between the x coordinate transmitted and 
reflected wave packets, and the y component wave packet. 

In Figs.[2|a-c), we illustrate solutions obtained using the reconstruction procedure 



(18 1, plotting the wavefunction on position space, \ip(x, y, t)| 2 , at three instants of 
time t — 0, 60, 100. In addition, to ease comparison, Fig. [2ji shows a plot that 
combines all three density distributions. The most obvious effect is the squeezing of 
the transmitted and the broadening of the reflected wave functions after the collision 
with the potential barrier. To understand this phenomenon we must realize that the 
energy band curvature induced by m(p y ) is associated to a momentum-dependent 
group velocity 

u>Px \Jc P x + m{p y yc l 

which in turn leads to a spreading of wavepackets in position space. This spreading 
is most relevant for small values of p x , that is when the particle crashes against the 
potential barrier. At this point, which is when the particle splits into a transmitted 
and a reflected component, a narrow band of p x components will be filtered and 
become part of the antiparticle branch, moving on with very little spreading and 
almost uniform group velocity. The reflected component, on the other hand, will get 
to see the whole curvature of the energy band and become very broad very quickly, 
until it reaches relativistic velocities (v g ~ c) and stops spreading. 

For performing the protocol, we envision to follow the lines developed for the 
1+1 case [12JI13] and use two trapped ions. One of them, together with two motional 
modes, will codify the spinor, including the x and y dependence of the wave packet. 
The other one will produce the linear potential in x. To be more specific, we consider 
two trapped ions trapped in a linear Paul trap, with two motional modes, with 
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frequencies v cm — v and v st — \fiv. To implement the Dirac-Klein dynamics given 
by Eq. (15), we proceed in an analogous way as in Refs. [9J [TU1 HH US] and use 
red and blue sidebands for each of the terms caiPi, i = x,y, as well as an AC-Stark 
shift for term mc 2 a z . The ax potential term will be implemented with a red and 
blue sidebands, with appropriate phases, applied to the second ion [T^ US]. The 



corresponding trapped-ion Hamiltonian will read 

r Kl,2+l 



D 



2n 



cm^cm^cm^.lPcm 



+ 2rj st A st n st a yA p st + Mia 



2,1 



(20) 



which coincides with Eq. (151 when one performs the analogies 2?7 cm A cm f2 cm = 
2?7 s tA s tf2st ^ c, Ml -H- mc 2 , and 77 cm Slo/A cm f-> a, when ion 2 is initialized in the 
eigenstate + of a x ^- 

For measuring the wave packet after performing the quantum simulation, one 
would proceed by mapping the motional state to the internal degrees of freedom 
of ion 2, and subsequently detecting the internal state by fluorescence detection, as 
exposed in Refs. [301 [SH 13 HD] - In these references, it was shown how to perform 
tomography of the ion motional state in an efficient way. 
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Figure 1. Probability \t/)(x,p y , t)\ 2 as a function of x, and the parameter p y 
(constant of the motion), for t = 0, 20, 40, 60, 80, 100 (ft = c = a = 1, m = 0.5). 

We point out that, when considering a nontrivial potential V(x, y) in x and y 
coordinates, say V{x,y) = a x x 2 + a y y 2 , which could be implemented with three ions 
in the dispersive limit, its simplification into a set of 1+1 problems cannot be done 
anymore, and the full bidimensional problem has to be considered. This means that, 
for reasonable sizes of motional Hilbert spaces of around 200 phonons per mode, one 
would have a Hilbert space dimension of 2 x 200 for the 1+1 dimensional case. On 
the other hand, for the 2+1 case considered here, one would have a Hilbert space 
dimension of 8 x 10 4 . The size of the states and the size of the associated Hamiltonian 
matrix force us to adopt particularly efficient schemes in the simulation, such as using 
finite differences or Fourier representations for the Hamiltonian and Trotter methods 
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Figure 2. Probability \ip(x, y)\ 2 as a function of x, y, for t = 0,60, 100, and a 
joint figure with the three times in a contour plot, (ti = c = a = 1, m = 0.5). 

for the unitary evolution. However, none of these techniques are without errors. They 
are severely affected by discretizations of time, and also by the boundary conditions, 
and even the simulations presented here are well within the limits of what can be 
accurately simulated, both in time and in space. These are the reasons why an actual 
experimental simulation of the 2+1 Klein scattering may have such an important 
value: without these constraints — memory, computational power, narrow boundary 
conditions — , it allows for an independent and perhaps more reliable verification of 
our predictions. 

4. Simulation of Dirac particles in magnetic fields 

So far, we have seen examples in which the trapped-ion simulation may be used to 
learn more about the Dirac equation, but the converse is also true: working in the 
quantum simulations may offer a new perspective on quantum optical problems. In 
particular, in this section, we will relate the simulation of Dirac particles in magnetic 
fields with the Jaynes-Cummings model, and use this to shed light on the topological 
features of the eigenstates of both models — the simulated one and the one used to 
implement the simulation. 

The topological properties of free Dirac fermions has become a very active research 
area in the last years [S3] . Roughly, in this context the energy band of a 2D free Dirac 
fcrmion can be seen as a mapping from momenta, {k x ,k y ) to a particular spinor 
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state on the sphere, S £ S, obtained from the Dirac wavefunction. The topologically 
nontrivial models are those for which this mapping winds up one or more times on the 
same sphere. The intriguing result is that such models, when embedded on surfaces 
with boundaries, present robust topologically protected "edge states", that is states 
that may transport charge or spin even when the bulk is placed in an insulating state 
(for instance with a large mass). 

We will see that some of these ideas may be cast in the context of trapped-ion 
simulations. For that, let us recall the model of a Dirac particle moving on a uniform 
two-dimensional magnetic field [M] , 

H = ca x p x + ca y [p y - eA(x, y)} + mc 2 a z . (21) 

We may choose an axial gauge in which A(x,y) = Bx. In this setup, the momentum 
along the y direction becomes a constant of motion, \p y , H] — 0. Changing units so that 
eB = 1, and displacing the "x" coordinate by an appropriate amount, x+p y — > x, we 
may rewrite the previous model as a simple Jaynes-Cummings model with a detuning, 

mc p x + ix \ fa, + f 2 



H = c\ . = cV2(a + a + a^a-) + mc z a z . (22) 

V Px — ix —mc J 

This model has a ladder of discrete energy eigenstates, 

E ±n = ±\/2c 2 n 2 + m 2 c A , (23) 

the so called Landau levels. The wavefunctions of these levels are confined along "x" 
and form plane waves along the "y" direction 

9 ±n (x,y) = ±( R a±n , M i, )e**y, n = l,2,... (24) 



^2 V P±n4>n-l(x) 

where (j) n (x) is the n-th level of a quantum harmonic oscillator and the coefficients 
(a± n , /3± n ) are the corresponding eigenvectors of the matrix 



*«=c( )■ (25) 



—iy/2n —mc 

At this point, instead of focusing as usual on the topology of energy bands, we 
will try to understand the topological properties of these individual states. For that 
we introduce the Wigner function of the spinor, given by 

W(x,p) = — / +s/2)* t (f + s/2)exp(-isp/h)d 2 s. (26) 

TTh J 

This quasiprobability distribution may be interpreted as a mapping from phase space 
(x,p) to points on the Bloch sphere, S £ S, 

s(x, Px ) = ~Tr{W(x,p x )a}, (27) 

where N — \\Tr{Wa}\\ is a normalization factor and we have dropped the y coordinate. 

Figure [3] shows the Wigner function of the two simplest Landau levels, n = ±1 and 
n = ±2, decomposed into two quantities, the polar angle (f> — tan -1 (s y /s x ), and the 
separation from the equator, s z . For the n — state we see that the pseudospin points 
North in the center of the phase space and rotates to the South at r — x 2 + p 2 — > oo. 
This is equivalent to a single covering of the sphere, or skyrmion, with a topological 
charge]!] v = 1. To achieve greater complexity we have to increase the number of 

% We do not take the orientation into account when computing this number, and instead define the 
charge as the absolute number of coverings. 
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Figure 3. Wigner function of the n = and n = 1 Landau levels, (a-b) 
and (c-d) respectively. The upper row shows the angle on the {S x ,S y } plane, 
tan - ((Sy)/(S X )), while the lower row plots the (S z ), along a P = cut. Note 
how the n = case implements a simple covering of the sphere, going from the 
north to the south poles, at x = and x = oo, while the n = 1 level covers the 
sphere twice, reaching the South pole at x = and x = oo. 



quanta, moving on to state n = 2, which performs a double covering of the sphere. 
This is done via an extended singularity of the angle <j) at radius r — 1, which is 
when the spin points North. In general we have verified that the number of coverings 
increases linearly with the number of quanta, v = ±n, leading to a family of topological 
defects with unprecedented richness. 

It may seem that these topological defects are mere artifacts of the Dirac 
or Jaynes-Cummings equations and that under normal circumstances, dephasing, 
spontaneous decay or other perturbations, will rapidly disappear. However, this is 
not the case. Pure dephasing, for instance along the "z" direction, is equivalent to a 
shrinking of the transverse components 

s{t) = {s x e-i\s y e-<\s z ). (28) 

This map contracts the total surface in Bloch space, which now looks closer to a rugby 
ball instead of a sphere. However, after suitable rescaling via a new normalization 
factor in Eq. (271, one finds that the total solid angle covered is the same and the 
charge, v, is not modified. 
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It is even more surprising to see that spontaneous emission has also a mild effect 
on the winding number. Applying the positive map that corresponds to this process 
onto the ion state, we see that the winding on the {S Xl S y } plane remains unaffected 
by the same reasons as before, and the only change is a shift of the weights of the 
different harmonic oscillator eigenstates, or equivalently s z . In Fig. [3ji, we summarize 
that process for the n = 2 Landau level, showing how dissipation tries to "unwrap" the 
sphere, removing one layer and falling down to a lower winding number. Note how, 
despite the effort in unwrapping the sphere, the defects remained pinned and the total 
number constant even for long times. One may consider a final form of losses: total 
dissipation in the form of friction in "x" or "p" . In either case, we expect population 
being transferred down the ladder of harmonic oscillator states, which this time would 
lead to a continuous unwrapping of the sphere v = n, n — 1, n — 2 . . . 0. 

Regarding the experimental implementation, the Hamiltonian for the Dirac 
particle in the field has, as we have shown, a direct reinterpretation in terms 
of a detuned Jaynes-Cummings model (22 1, which makes its simulation rather 
straightforward. Measurement of the Wigner function is also possible using standard 
trapped ion techniques [49] , but now the reconstruction has to be performed without 
tracing out the internal state of the ion. To do this will require using an ancillary 
ion for measurement and a different ion to implement the spinor. The ancillary ion 
will remain in the ground state until the end of the experiment, when it will be 
used to gather statistics about the vibrational mode that implements the Jaynes- 
Cummings Hamiltonian [49]. The novelty is that this statistics will have to be 
performed conditionally on the spinor-ion, which should be simultaneously measured 
in one of the a x , a y or a z basis, as it is usual for wavefunction reconstruction. The 
resulting statistics and accuracy should be more than enough to allow for accurate 
reconstruction of the v = 1 and v = 2 skyrmions. 



5. Two Dirac particles interacting via a classical potential: Simulation of 
an analogue of MIT bag model 

The quantum field theory for the strong interaction, known as quantum 
chromodynamics (QGD), accurately predicts the behavior of nuclei, hadrons, and 
quarks and gluons. On the other hand, for low energies and long distances the 
theory becomes non-perturbative (because the coupling becomes too large) such that 
Feynman-diagram expansion cannot be applied. 

Many decades ago, some effective theories were developed that qualitatively 
and quantitatively described the physics of the strong interaction even in the non- 
perturbative regime. Among those theories, arguably, the best known is the MIT 
bag model [46] |47] . In its simplest version, it consists of two or more Dirac particles 
contained inside a volume, a bag, whose energy grows linearly with the volume. In 
this way, this model can simulate behavior like asymptotic freedom and confinement: 
when the fcrmions go far away from each other, the energy grows steadily and there 
is a tendency to come back to the original position, i.e., confinement takes place. On 
the other hand, when the fermions are nearby, the potential energy becomes negligible 
and the fermions behave as if they were free, i.e., asymptotic freedom holds. Although 
these models can be solved in a classical computer, and are phenomenological, they 
have quite a good agreement with experiments and predictive power. Even nowadays, 
this kind of methods, either with color charge or without it, string models, etc., are 
the basic tools to analyze QCD at low energies in the nonperturbative regime, with 
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semianalytic tools (without resorting to lattice gauge theory, which requires huge 
computational power) [55] . We believe it would be interesting to implement analogues 
of the MIT bag model in a quantum simulator, even if it is just for reproducing 
the physics of quark confinement and asymptotic freedom in a table top tunable 
experiment. In addition, quantum simulations of these models with state-of-the- 
art technology could already go beyond what can be computed classically (for large 
motional Hilbert space dimensions). In this Section, we propose the simulation of 
a simplified analogue of MIT bag model in 1+1 dimensions. Thus, we propose a 
simulation of two Dirac equations, one for each Dirac particle, that are coupled by 
a potential which grows with the distance, in our case quadratically for the sake of 
experimental feasibility. This model could describe a quark and antiquark coupled 
by gluons inside a meson, and cannot be solved analytically. By appropriately 
tuning the lasers in an experiment, one could perform phase transitions between the 
asymptotically free and the confined meson phases, like the one that supposedly took 
place in the early Universe. 

The model we aim to simulate is 

Habm = ca Xj iPi + ca x ^p 2 + mc 2 a yA + mc 2 a y ^ + V (xi - x 2 ) 2 . (29) 

In order to simulate two Dirac particles (e.g., quark-antiquark constituting a meson), 
each obeying a 1+1 Dirac equation and with a potential that grows with the separation 
between the particles, we envision to use three ions: two of them will codify the positive 
and negative energy spinors of the two Dirac particles, while the third one will be used 
to generate the attractive potential among them. The normal modes for a three ion 
system with equal mass and charge are 

Q cm = -^(xi+x 2 + x 3 ), P cm = ~^=(Pi +P2 +P3), (30) 

Qr = - -y=(Xx-X 3 ), Pr = -^(Pl-P3), (31) 

Q 3 = ^(zi -2x 2 + x 3 ), P 3 = -±={ Pl -2p 2 +p 3 ). (32) 

We will consider the two first modes to codify the two free Dirac equations plus the 
potential. By considering the Hamiltonian 

-ffABM.sim = 2r/ cm A cm fl cm (<J x _ 1 — <J x ^)P cm + 2"q r A r tt r (<T Xv i + <J x ^)Pr 

+ ft£la yi i + HQa Vj3 + fejggi— - a Xi2 Qem + ^A 3 er 2i2 , (33) 

we simulate an equivalent dynamics to two free Dirac equations with a potential in 
the relative coordinate, 

H abm = cct x aPi + ca Xi3 p 3 + mc 2 o- yA + mc 2 cry t3 + V (xi - x 3 ) 2 (34) 

P 

= c(a Xi x - a Xt3 )p r + c(a Xjl + (r Xl a)-~ + mc 2 a yA + mc 2 a Vi3 + V x 2 , 

where p r = (jp\ —p 3 )/2, and P cm = p\ +p 3 are the relative and center of mass momenta 
of the system of two particles that we want to simulate, and x r — x\ — x 3 is the 
relative coordinate. To make the analogy complete, one has to make the substitutions 

Pcm Pr: Pr ^ Pan/2^ Q cm ^ X ri with t 2f] c7n A cnl d cm 2Tj r /S. r Cl r C, H£l TTIC , 

and the potential is obtained in the dispersive limit, for hA 3 ^> hrj cm Vl, 3 (a\ m + a cm ), 
considering ion 2 in the + eigenstate of er Z:2 , and with Vq = (hn cm Q 3 / A cm ) 2 / hA 3 . 
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With this simulation, we may analyze confinement and asymptotic freedom 
by playing with Vq potential. One would expect in principle, though predictions 
are complicated in systems with no analytical solutions, to have two independent 
dynamics. One of them will happen for the case in which the fermions are nearby 
(x r ~ 0) and the momentum \p r \ is small (asymptotic freedom). And another more 
complicated dynamics, coupling the two fermionic motions, for those cases in which the 
two fermions fly away (x r large) (confinement). Additionally, when Vq(x^) > 2mc 2 , 
one would expect Klein tunneling to take place. In this model, we interpret this Klein 
tunneling as a pair creation of a new quark and a new antiquark: when the original 
quark and antiquark go very far away from each other, there is energy available for 
creating a new antiquark (that gets bound to the original quark) and a new quark 
(that gets bound to the original antiquark), such that the original pair can split up 
without violating the color neutrality standard in QCD for free composite particles. 
This is the natural interpretation, given that the usual explanation for Klein tunneling 
in quantum field theory is the creation of a particle/antiparticle pair. In QCD, the 



Relativistic quantum mechanics with trapped ions 



15 



splitting of a hadron into new hadrons is well known and takes place, for example, in 
the jet emission in high-energy colliders. 

We have realized numerical simulations of this model, and we show our results 
in Fig. [4] There, we plot the probability \ifr(x r , t )\ 2 as a function of x r and t, for 
the two Dirac fermion system as given by Eq. (34 1, for P cm = 2 (a constant of 
motion), mc 2 = 1, and Vq — 0.5, for the initial values (a) (p r ) = 0, II = +1, (b) 
(p r ) = 2, II = +1, (c) (p r ) = 0, n = -1, and (d) (p r ) = 2, II = -1. In this notation, 
II denotes the expected value of the operator ft = \{(J X ,\ ~ ^1,3) on the initial state. 
We consider initial Gaussian states in x r , normalized in this coordinate, and constant 
P cm , that will remain unchanged during the evolution, for the sake of simplicity and 
ease of the numerical simulation. 

Despite the fact that this problem is not analytically solvable, and it is difficult to 
get a clear intuitive prediction of the behavior that Eq. (34 1 will produce, we observe 
interesting features. For II = +1, and (p r ) = the dynamics is quasifree, besides a 
small Klein tunneling appearance. This may be a signature of asymptotic freedom, 
because the Dirac fermions do not couple much, as seen in (a). For II = +1 and 
(Pr) = 2, the dynamics of the two-fermion system is more involved: the wave function 
evolves in a more complex way for larger relative momentum, splitting up, bouncing 
back and getting distorted. This may be related to confinement taking place, when 
the two fermions try to separate, as seen in (b). On the other hand, for II = — 1 
(another spinor initial state), we observe that, for finite relative momentum, (p r ) = 2 
(d), the two fermions have a larger Klein tunneling than for zero relative momentum, 
(p r ) = (c). This may also indicate that, for large enough relative momentum, the 
Dirac fermion pair has enough energy to create another quark-antiquark pair and 
escape from each other while keeping color neutrality (d). For relative momentum 
equal to zero, there is little Klein tunneling (c) associated to asymptotic freedom. 

Summing up, we have shown than an experiment with three ions can simulate 
an analogue of the well-known MIT bag model [46l [47] for two Dirac particles in 
1+1 dimensions. Despite their simplicity, such experiments would allow us to study 
signatures of asymptotic freedom and confinement. In addition, those experiments 
would easily reach simulation regimes which are beyond current numerical simulation 
possibilities. For example, here, we were severely limited by the problem size 
and complexity, preventing us from studying things like adiabatic preparation of 
eigenstates, a full diagonalization of the spectrum, or longer time dynamics. 



6. Conclusions 



In this article, we have shown that trapped-ion physics provides a flexible platform to 
simulate a number of interesting cases and effects in relativistic quantum mechanics. 
In particular, we have studied bidimensional relativistic scattering, topological effects 
of Dirac particles in magnetic fields, and two Dirac particles coupled by a potential, 
an analogue of nuclear physics bag models. 

Summarizing, we have described the bidimensional relativistic scattering for an 
^-dependent potential, showing that this case can be solved with the help of the 
one-dimensional case. In addition, we showed it contains interesting novel features, 
like entanglement between transmitted/reflected wave packets and the transverse 
momentum. We have also given a heuristic analysis of the complexity of more 
complicated potentials that will already approach the classical computational limit 
of resources. We have introduced a relativistic quantum mechanics model -a Dirac 
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particle in a magnetic field- that could shed light and establish interesting analogies 
with an equivalent quantum optical model, namely, the Jaynes-Cummings model 
with detuning. We showed a possible quantum simulation of this case, including 
its topological features in trapped ions. Finally, we have shown that three trapped 
ions together with two motional modes suffice to implement the dynamics of two 
Dirac particles that are coupled by a confining potential, like in the bag models of 
nuclear physics. This could allow the implementation in a table top experiment of 
quark confinement and asymptotic freedom associated to QCD, while moving from 
one regime to the other by switching on and off a laser. This would amount to a 
phcnomenological simulation of the phase transition between the quark-gluon plasma 
and the hadronic confined phases that took place in the early universe. 

In a physical implementation, strings of Ca + ions, for instance, could be used 
and manipulated with laser light. The spinor degrees of freedom can be encoded in 
long-lived electronic states of this ion, as was done before in quantum simulations 
of relativistic quantum mechanics [101 113] . Such states can have coherence times of 
several ms, while coupling strengths of a few 100 kHz between the spinor states can 
be obtained. For each motional degree of freedom, a normal motional mode of the ion 
string can be used and manipulated with lasers coupling the internal and motional 
degrees of freedom. The implementation of the proposed quantum simulations requires 
the validity of the Lamb-Dicke approximation, which for Ca + ions in a trap with 
secular frequency of ~ 1 MHz and manipulated with 729 nm laser light is well justified 
(Lamb-Dicke parameter r\ ~ 0.05 -C 1) [33] . In recent experiments [13], it has been 
shown that the motional state of the ions can be coherently manipulated while reaching 
states that have >100 phonons on average, and the coherence can be sustained for 
more than 10 ms. These parameters could be further improved by using, e.g., ions 
that are less sensitive to decoherence sources and by further decreasing the Lamb-Dicke 
parameter with the choice of a different beam geometry. 

We believe that quantum simulation is one of the most promising fields inside 
quantum information science. At the same time, among other research lines, quantum 
simulations of relativistic quantum mechanics may soon allow us to consider problems 
that are difficult or even impossible for classical computers. In this sense, quantum 
simulators will give us the opportunity to enter into unexplored regimes of the physical 
world. Trapped ions offer one of the most promising platforms for achieving this goal. 
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